function viewSizes( sid, part, ebsd, saveres, thrd, mgs)
%View information  about grain structure of EBSD.
%   Draw orientation maps
%
% Syntax
%   viewSizes( sid, part, ebsd, saveres, thrd, mgs)
%
% Input
%   sid     - sample id: 's01', 's02', 's03', 's04', 't01', 'p01', 'p02'
%   ebsd    - EBSD data if 0, try load useing function "[sid '_load']"
%   saveres - save image to disk
%   thrd    - boundary detection threshold in degree
%   mgs     - grain's size threshold (lesser grains deleted)
%
% History
% 07.12.12  Add description of the function. Rename 'ma' to 'thr', 's' to
%           'mgs'.
% 14.04.13  Add saveing of comment.
% 16.04.13  Separate degree and radian 'thr'.
% 16.05.13  Add automatic intersect method.

% Convert to radian
thr = thrd*degree;

OutDir = checkDir(sid, 'sizes', saveres);

ebsd = checkEBSD(sid, ebsd, 'Fe');

comment = getComment();

% File name prefix
prefixg = [sid '_' num2str(part) '_grains_' num2str(thrd) '_' num2str(mgs) '_x'];
prefixh = [sid '_' num2str(part) '_hist_' num2str(thrd) '_' num2str(mgs) '_x'];

% grains = getGrains(ebsd, thr, mgs);
% 
% figure();
% plot(grains,'antipodal',1);
% saveimg( saveres, 1, OutDir, prefixg, '', 'png', comment );
% 
% % Grains' area
% ar = area(grains);
% 
% % Grains' equivalent radius
% er = equivalentradius(grains);
% 
% % Grains' shape factor
% sf = shapefactor(grains);
% 
% % Grains' aspect ratio
% as = aspectratio(grains);
% 
% disp(['Grains size information for ' prefixg ]);
% grainStat(grains, saveres, OutDir, prefixh);
% 
% % Logarithm scale for grains' area
% bins = exp(-1:0.4:log(max(ar)));
% 
% % Plot area histogram
% figure();
% [nlog, xlog] = hist(ar,bins);
% bar( xlog, nlog );
% saveimg( saveres, 1, OutDir, prefixh, 'area', 'png', comment );
% 
% savexy( xlog, nlog, saveres, OutDir, prefixh, 'area_log', 'Area Disribution (Log)', comment );
% 
% [y, x] = hist(ar);
% savexy( x, y, saveres, OutDir, prefixh, 'area', 'Area Disribution', comment );
% 
% savexy( 1:length(ar), ar, saveres, OutDir, prefixh, 'area_raw', 'Area RAW', comment );
% 
% 
% % Plot equivalent radius histogram
% figure();
% bins = [0.25 0.5 1 2 4 8 16 32];
% bar( hist(er,bins) )
% set(gca,'XTickLabel',{'0.25';'0.5';'1';'2';'4';'8';'16';'32'})
% title('Radius');
% saveimg( saveres, 1, OutDir, prefixh, 'radius_log', 'png', comment );
% 
% % Plot equivalent radius histogram
% figure();
% bins = [0.25 0.5 1 2 4 8 16 32];
% bins = pi*(bins.*bins);
% bar( hist(ar,bins) )
% set(gca,'XTickLabel',{'0.25';'0.5';'1';'2';'4';'8';'16';'32'})
% title('Area');
% saveimg( saveres, 1, OutDir, prefixh, 'radius_log_ar', 'png', comment );
% 
% % Plot equivalent radius histogram
% figure();
% bins = [2:2:20];
% bar( bins, hist(er,bins) )
% saveimg( saveres, 1, OutDir, prefixh, 'radius_lin', 'png', comment );

%% Grain size vs threshold angle
% Output diectory for temporary grains map
OutDirSec = checkDir(sid, 'sizes\sec', saveres);

% Average grain size for different thr
alist = 2:5:60;
nn  = length(alist);
ags = zeros(1,nn);
aga = zeros(1,nn);
agaa = zeros(1,nn);
agr = zeros(1,nn);

if  any(cellfun(@(x) strcmpi(x,'ci'),get(ebsd)))
    ebsd = ebsd(get(ebsd,'CI') > 0.1);
end

% Grain size from boundary length
grains = getGrains(ebsd, thr, mgs);
mori = calcBoundaryMisorientation(grains);
ang = angle(mori)/degree;

n = zeros(1,nn);
for i = 1:nn
    n(i) = sum(ang > alist(i));
end
[dx,dy,~,~,sx,sy] = getStep(ebsd);
dd = (dx+dy)/2;
n = n*dd/1.273;
agl = 2*(sx*sy)./n;

% Grain size for different thr angle
for i = 1:nn
    a = alist(i);
    try
        [ags(i), aga(i), agaa(i), agr(i)] = difGrainSize(ebsd, a, mgs, OutDirSec, prefixg);
    catch err
        warning(err.message);
    end
end

figure;
plot(alist, 1.27*ags, alist, 2*sqrt(aga/pi), alist, 2*agr, alist, agaa, alist, agl);
title('Grain size vs threshold angle');
legend('intersect "1.27*int"', 'from area "2*sqrt(area/pi)"', 'from eq radius "2*er"', 'from areafraction', 'from misori', 'Location','NorthWest');
saveimg( saveres, 1, OutDir, prefixg, 'thr', 'png', comment );

% Save data 
saveopt( alist, [1.27*ags; 2*sqrt(aga/pi); 2*agr; agaa; agl]',...
    {'intersect', 'from area', 'from eq radius', 'from areafraction', 'from misori'},...
    saveres, OutDir, prefixg, 'thr', comment );
end


function [ags, aga, agaa, agr] = difGrainSize(ebsd, thrd, mgs, OutDir, prefix, saveres)

    comment = getComment();

    % Calculate grain map
    grains = getGrains(ebsd, thrd*degree, mgs);
%     grains = smooth(grains,1);
    
    % Plot map and save it
    figure;
    plotBoundary(grains); axis off;
    saveimg( 0, 1, OutDir, prefix, ['sec' num2str(thrd)], 'png', comment );
    close all;
    
    % Generate grain map image filename
    fname =  [OutDir '\' prefix '_sec' num2str(thrd) '.png'];
    
    % Get ebsd data size for calculate image resolation (um per pixel)
    [~,~,~,~,sx,sy] = getStep(ebsd);
    
    % Run intersect method
    ags = grSizeSec( fname, 150, [sx sy], 'crop', 'grid' );
    
    % Grains' area
    ar = area(grains);
    dm = ceil(2*sqrt(max(ar)/pi));
    dd = 0:dm;
    aga = mean(ar);
    af = areafraction(ar, dd);
    agaa = sum(af.*dd)/sum(af);
    
    % Grains diameter
    r = sqrt(ar/pi);
    agr = mean(r);
    
%     % Debug
%     figure;
%     em = ceil(log(max(ar)));
%     edges = exp(0:1:em);
%     h = histc(ar,edges);
%     bar(h);
%     
%     figure;
%     em = ceil(log(max(r)));
%     edges = exp(0:1:em);
%     h = histc(r,edges);
%     bar(h);
%     
%     figure;
%     hist(r);
%     
%     figure;
%     bar(dd,af/sum(af)*100,'histc')
end
